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Abstract Algorithms for the symbolic computation of polynomial conservation 
laws, generalized symmetries, and recursion operators for systems of nonlinear 
differential-difference equations (DDEs) are presented. The algorithms can be used 
to test the complete integrability of nonlinear DDEs. The ubiquitous Toda lattice il- 
lustrates the steps of the algorithms, which have been implemented in Mathematica. 
The codes InvariantsSymmetries.m and DDERecursionOperator.m can 
aid researchers interested in properties of nonlinear DDEs. 
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1 Introduction 

A large number of physically important nonlinear models are completely integrable, 
i.e., they can be linearized via an explicit transformation or can be solved with the 
Inverse Scattering Transform. Completely integrable continuous and discrete mod- 
els arise in many branches of the applied sciences and engineering, including clas- 
sical, quantum, and plasma physics, optics, electrical circuits, to name a few. Math- 
ematically, nonlinear models can be represented by ordinary and partial differential 
equations (ODEs and PDEs), differential-difference equations (DDEs), or ordinary 
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and partial difference equations (OziEs and PziEs). This paper deals with integrable 
nonlinear DDEs. 

Completely integrable equations have nice analytic and geometric properties re- 
flecting their rich mathematical structure. For instance, completely integrable PDEs 
and DDEs possess infinitely many conserved quantities and generalized (higher- 
order) symmetries of successive orders. The existence of an infinite set of gener- 
alized symmetries can be established by explicitly constructing recursion operators 
which connect such symmetries. Finding generalized symmetries and recursion op- 
erators is a nontrivial task, in particular, if attempted by hand. For example, in Q 
and lfl4l an algorithm is presented to compute recursion operators for completely 
integrable PDEs, which was only recently implemented in Mathematica (TJ. 

Based on our earlier work in ff7], []9], and iTTOl . we present in this paper algo- 
rithms for the symbolic computation of conserved densities, generalized symme- 
tries, and recursion operators of nonlinear systems of DDEs . Such systems must be 
polynomial and of evolution type, i.e., the DDEs must be of first order in (continu- 
ous) time. The number of equations in the system, degree of nonlinearity, and order 
(shift levels) are arbitrary. Furthermore, the current algorithms only cover polyno- 
mial densities, symmetries, and recursion operators. 

We use the dilation (scaling) invariance of the system of DDEs to determine 
the candidate density, symmetry, or recursion operator. Indeed, these candidates are 
linear combinations with undetermined coefficients of scaling invariant terms. Upon 
substitution of the candidates into the corresponding defining equations, one has to 
solve a linear system for the undetermined coefficients. After doing so, the coeffi- 
cients are substituted into the density, symmetry, or recursion operator. If so desired, 
the results can be tested one more time, in particular, by applying the recursion op- 
erators to generate the successive symmetries. 

If the system of DDEs contains constant parameters, the eliminant of the lin- 
ear system for the undetermined coefficients gives the necessary conditions for the 
parameters, so that the given DDEs admit the required density or symmetry. In anal- 
ogy with the PDE case in |8|], the algorithms can thus be used to classify DDEs with 
parameters according to their complete integrability as illustrated in (9) and iflOl . 

As shown in [4j, once the generalized symmetries are known, it is often possible 
to find the recursion operator by inspection. If the recursion operator is hereditary, 
as defined in |6), then the equation will possess infinitely many symmetries. If, in 
addition, the recursion operator is factorizable then the equation has infinitely many 
conserved quantities. 

Computer algebra systems can greatly help with the search for conservation laws, 
symmetries, and recursion operators. The algorithms in this paper have been imple- 
mented in Mathematica. The computer codes (see fPH l. can be used to test the 
complete integrability of systems of nonlinear DDEs, provided they are polynomial 
and of first order (or can be written in that form after a suitable transformation). 

With InvariantsSymmetries.m, in Q, J9), and Q0|, Goktas and Hereman 
computed polynomial conserved densities and generalized symmetries of many 
well-known systems of DDEs, including various Volterra and Toda lattices as well 
as the Ablowitz-Ladik lattice (for additional results and references, see, e.g., |[T6l ). 
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The existence of, say, a half dozen conserved densities or generalized symmetries 
is a predictor for complete integrability. Finding a recursion operator then becomes 
within reach. An existence proof (showing that there are indeed infinitely many 
densities or generalized symmetries) must be done analytically, e.g., by explicitly 
constructing the recursion operator which allows one to generate the generalized 
symmetries order by order. Numerous explicit examples have been reported in the 
literature but novices could start with the book by Olver lTT8l to learn about recursion 
operators for PDEs. To alleviate the burden of trying to find a recursion operator by 
trial and error, we present a new Mathematica program, DDERecursionOper- 
ator.m, based on the algorithm in Section 5. Like InvariantsSymmetries.m, 
after thorough testing, DDERecursIOnOperator.M will be available from |[T2l . 

If one cannot find a sufficient large number of densities or symmetries (let alone, 
a recursion operator), then it is unlikely that the DDE system is completely inte- 
grable, at least in that coordinate representation. However, our software does not 
allow one to conclude that a DDE is not completely integrable merely based on the 
fact that polynomial conserved densities and generalized symmetries could not be 
found. Polynomial DDEs that lack the latter may accidentally have non-polynomial 
densities or symmetries, or a complicated recursion operator, which is outside the 
scope of the algorithm described in Section 5. 

Currently, our algorithm fails to find recursion operators for the Belov-Chaltikian 
lattices El EQ] EE] and lattices due to Blaszak and Marciniak 131 l20l |2T1 |24l . In the 
near future we plan to generalize the recursion operator algorithm so that it can 
cover a broader class of nonlinear DDEs. 

The paper is organized as follows. Basic definitions are given in Section 2. In 
Section 3, we show the algorithm for conservation laws, using the Toda lattice as 
an example. Using the same example, Sections 4 and 5 cover the algorithms for 
generalized symmetries and recursion operators, respectively. In Section 6, we draw 
some conclusions and briefly discuss future research. 



2 Key Definitions 

Consider a system of nonlinear DDEs of first order, 

U;i = F(u„_{<,...,U„_i,U„,U„ + l,...,U„ + m), (1) 

where u„ and F are vector-valued functions with N components. This paper only 
covers DDEs with one discrete variable, denoted by integer n, which often corre- 
sponds to the discretization of a space variable. The dot stands for differentiation 
with respect to the continuous variable (often time ?). Each component of F is as- 
sumed to be a polynomial with constant coefficients. If parameters are present in 
(Q]l, they will be denoted by lower-case Greek letters. F depends on u„ and a finite 
number of forward and backward shifts of u„. We denote by £ (m, respectively), the 
furthest negative (positive, respectively) shift of any variable in the system. Restric- 
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tions are neither imposed on the degree of nonlinearity of F, nor on the integers I 
and to, which measure the degree of non-locality in (|T). 



2.1 Leading Example: The Toda Lattice 

One of the earliest and most famous examples of completely integrable DDEs is the 
Toda lattice, discussed in, for instance, fl22l : 

% = exp (y„_ i -y n )- exp (y n -y n+ \), (2) 

where y n is the displacement from equilibrium of the «th particle with unit mass 
under an exponential decaying interaction force between nearest neighbors. In new 
variables (u n , v„), defined by u„ =j„, v„ = exp(y„ — y„+i ), lattice (fSJ) can be written 
in polynomial form 

v„ = v n (u n -u n+ i). (3) 

The Toda lattice (O will be used to illustrate the various algorithms presented in 
subsequent sections of this paper. 



2.2 Dilation Invariance 

A DDE is dilation invariant if it is invariant under a dilation (scaling) symmetry. 
2.2.1 Example 

Lattice (0 is invariant under scaling symmetry 

(f,M„,v„) ->• (A _1 f,A 1 M„,A 2 v„). (4) 



2.3 Uniformity in Rank 

We define the weight, w, of a variable as the exponent of the scaling parameter (A) 
which multiplies that variable. Since A can be selected at will, t will always be 
replaced by t and, thus, w(^) = w(D,) = 1. 

Weights of dependent variables are nonnegative, rational, and independent of n. 
For example, w(u n -i) = ■■■ = w(u n ) = ■ ■ ■ = w(u n+ 2}- 
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The rank, denoted by of a monomial is defined as the total weight of the 
monomial. An expression is uniform in rank if all of its terms have the same rank. 

Dilation symmetries, which are special Lie-point symmetries, are common to 
many DDEs. Polynomial DDEs that do not admit a dilation symmetry can be made 
scaling invariant by extending the set of dependent variables with auxiliary param- 
eters with appropriate scales as discussed in J5] and ifTOl . 



2.3.1 Example 

In view of (0}, we have w(u„) = 1, and w(v„) = 2 for the Toda lattice. In the first 
equation of (f3]), all the monomials have rank 2; in the second equation all the mono- 
mials have rank 3. Conversely, requiring uniformity in rank for each equation in (O 
allows one to compute the weights of the dependent variables (and, thus, the scaling 
symmetry) with simple linear algebra. Balancing the weights of the various terms 
of each equation in (01 yields 

w(m„) + 1 = w(v n ), 

w(v„) + l = w(u„)+w(v„). (5) 

Hence, 

w(u n ) = 1, w(v„) = 2, (6) 

which confirms (0). 



2.4 Up- Shift and Down-Shift Operator 

We define the shift operator D by Du„ = u„+i . The operator D is often called the up- 
shift operator or forward- or right-shift operator. The inverse, D~ 1 , is the down-shift 
operator or backward- or left-shift operator, D _1 u„ = u„_i. Shift operators apply 
to functions by their action on the arguments of the functions. For example, 

DF(u 

= F(Du„_(>, • • • , Du„_i , Du„, Du„+i , . . . , Du„ +m ) 

= F(u„_m-i,...,u„,u„ + i,u„ + 2,--- ,u„ +m+ i). (7) 



2.5 Conservation Law 



A conservation law of (Q3, 

D t p+AJ = 0, 



(8) 
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connects a conserved density p to an associated flux J, where both are scalar func- 
tions depending on u„ and its shifts. In ©, which must holds on solutions of ([TJ, D, 
is the total derivative with respect to time, A = D — I is the forward difference op- 
erator, and I is the identity operator. For readability (in particular, in the examples), 
the components of u„ will be denoted by u„ , v n , w„ , etc. In what follows we consider 
only autonomous functions, i.e., F,p, and J do not explicitly depend on t and «. 

A density is trivial if there exists a function y so that p = Ay. We say that two 
densities, pW and p^ 2 \ are equivalent if and only if p' 1 ' +kp^ = A\\f, for some 
y/ and some non-zero scalar k. It is paramount that the density is free of equivalent 
terms for if such terms were present, they could be moved into the flux J. 

Compositions of D or D _1 define an equivalence relation (=) on monomial 
terms. Simply stated, all shifted terms are equivalent, e.g., k„_iV„+i = u n v n+ 2 = 
u n +2V n +4 = «n-3V„-i since 

u n -iv n+ \ = u„v n+ 2 — 4(h„_iv„ + i) 

= U n+ 2V n+ 4 - A (m„+1 V„ + 3 + U n V n+ 2 + U n -lV n+ \) 

= U n -3V n -l + A(u„_ 2 Vn + U„-3V„-l)- (9) 

This equivalence relation also holds for any function of the dependent variables, 
but for the construction of conserved densities we will apply it only to monomial 
terms (/,•) in the same density, thereby achieving high computational efficiency. In 
the algorithm used in Section 3, we will use the following equivalence criterion: 
two monomial terms, t\ and ?2, are equivalent, t\ = ?2, if and only if t\ = D r t2 for 
some integer r. If t\ = ?2 then t\ = t% + AJ for some J dependent on u„ and its 
shifts. For example, u n -2U n = u n -\u n+ \ because M„_2"n = D~ u n -\u n +\. Hence, 

U„-2Un = U„-\U n+ \ + [—U„-\U n +\ + U n -2U n ] = u n -\u n+ y +AJ with / = —U n -2U n . 

For efficiency, we need a criterion to choose a unique representative from each 
equivalence class. There are a number of ways to do this. We define the canonical 
representative as that member that has (i) no negative shifts and (ii) a non-trivial 
dependence on the local (that is, zero-shifted) variable. For example, m„m„+2 is the 
canonical representative of the class 

{• • • j U n -2U n ,U n -\ M„+l , U n U n +2, U n+ 1 K„+3 , •••}• 

In the case of, e.g., two variables [u„ and v„), u n+ 2V„ is the canonical representative 
of the class 

{• • • ,U„-\ V„_3 , U n V n -2 , U n+ 1 V„_ 1 , U n+ 2V„ , U n+ 3V n+ 1 , • • • }. 

Alternatively, one could choose a variable ordering and then choose the mem- 
ber that depends on the zero-shifted variable of lowest lexicographical order. The 
code in lfl2l uses lexicographical ordering of the variables, i.e., u„ -< v n ~< w n , 
etc. Thus, m„v„_2 (instead of u n+ 2V n ) is chosen as the canonical representative of 

{■ ■ ■ ,U n -lV n -3,U n V n -2,U n+ lV n -l,U n+ 2V n ,U n+ 3V n+ l,- ■ ■ }. 
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It was shown in [ 17] that if p is a density then D k p is also a density. Hence, using 
an appropriate "up-shift" all negative shifts in a density can be removed. Without 
loss of generality, we thus assume that a density that depends on q shifts has canon- 
ical form p(u„,u„+i,- • • ,U n+q ). 



2.5.1 Example 

Lattice (0) has infinitely many conservation laws (see, e.g., ifTTl ). Here we list the 
densities of rank R < 4 : 

P [l) =u n , (10) 
P W = \u 2 n + Vn , (11) 
p (3) = ±i4 + K n (v„_i+v„), (12) 

P (4) = 3«n + "n(v„-l+V„) + M„M„ + iV n + iv^ + V„V„ + i. (13) 

The first two density-flux pairs are easily computed by hand, and so is 

pi 0) =ln(v„), (14) 
which is the only non-polynomial density (of rank 0) . 



2.6 Generalized Symmetry 

A vector function G(u n ) is called a generalized symmetry of (HJ if the infinitesimal 
transformation u„ — > u„ + eG leaves (fl]i invariant up to order e. As shown by fTSl . 
G must then satisfy 

D,G = F'(u„)[G] (15) 

on solutions of (Q]), where F'(u„)[G] is the Frechet derivative of F in the direction 
ofG. 

For the scalar case (N = 1), the Frechet derivative is 

d d F 

F'(u n )[G] = —F(u n + eG)| £=0 = £ ^ D k G, (16) 

de Y du n+k 

which, in turn, defines the Frechet derivative operator 




In the vector case with, say, components u„ and v n , the Frechet derivative operator 
is a matrix operator: 



F'(u„) 



Unal G6kta§ and Willy Hereman 



i v dF 2 r^k v dF 2 nk 



(18) 



Applied to G = (Gi G2) , where T is transpose, one obtains 

f ;'(u, ! )[G]=£^D*G 1 +£^-D*G 2 , (19) 

V du n+k V dv n+k 

with i = 1,2. In ( [ToT l and ( fT9l summation is over all positive and negative shifts 
(including k = 0). The generalization of ( fT8l to a component system is straight- 
forward. 



2.6.1 Example 

As computed in [Oj, the first two non-trivial symmetries of (O are 

Vn-V„-l 
V n (u n+ i — U n ) 

' v n {u n + u n+ \) - v„_i(« n _i + u n ) ' 
v«(w^ +1 - u\ + v n+ i - v n _i) 



G« = 



G (2) 



2. 7 Recursion Operator 

A recursion operator S? connects symmetries 

G u+s) =MG {] \ (20) 

where j = 1 , 2, • • ■ , and s is the gap length. The symmetries are linked consecutively 
if s = 1. This happens in most (but not all) cases. For A^-component systems, Si is 
anN xN matrix operator. 

With reference to [18) and fl23l . the defining equation for ffl is 

D t Si +[<%,¥' (u„)] 

= ^ + SH [F] + St o F' (o„ ) - F' (u„) o & = 0, (21) 
at 

where [ , ] denotes the commutator and o the composition of operators. The operator 
F'(u„) was defined in ( fT8l . S?'[F] is the Frechet derivative of 0t in the direction of 
F. For the scalar case, the operator 0£ is often of the form 
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^ = C/( Mn )^((D-I)- 1 ,D- 1 ,I,D) V(u n ), (22) 

and then 

M'[F] =Y,(D k F) 1 ^-0V + Y,U0{D k F)J^. (23) 

For the vector case, the elements of the N xN operator matrix M are often of the 
form 

Mij = Uij(u n ) ffij ((D - rrSD-Sl.D) Vy(u„). (24) 
Hence, for the 2-component case 

ae [f] y = E (d** ) en vn + L (d*f 2 ) i^L ^ v< . 

+ £ C/y tfy (D*Fi) ^ +L f/y <?y (D^ 2 ) J^L. (25) 

2.7.1 Example 

The recursion operator of (0 is 

M„I D 1 +1+ (v n - v„_i)(D— I) -1 7-I \ 

(26) 

v„I + v„D u n+ il + v n (u n+ i -w„)(D-l) _1 ^I / 
It is straightforward to verify that MG W = G (2) with G {1) in (|20j and G (2) in ©. 



3 Algorithm for Conservation Laws 

As an example, we will compute the density p' 3 ' (of rank R = 3) given in ( TT2b . 



3.7 Construct the Form of the Density 

Start from Y = {u n ,v n }, the set of dependent variables with weights. List all mono- 
mials in u and v of rank R = 3 or less: j$ = {u%,u%,u n v n ,u n ,v n }. Next, for each 
monomial in J£, introduce the correct number of f -derivatives so that each term has 
rank 3. Using ([3J, compute 
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dV 

— = U 3 
dt° 



n — ll n V ni 

dr 

— =2u n u„ =2u„v n -i -2u n v„, (27) 
dt 

dv„ 

— — = V„ = U n V n — U„+lV n , 

dt 

d 2 « n _ &Un _ d(v n -l - V' n ) 

dt 2 ~~ dt dt 

= u n -\v n -\ —u n v n _\ —u n v n + u n +iv„. 

Gather the terms in the right hand sides in (l27t to get 

y = {ul,u n v n -\, u n v n , u„_ i v„_ i , u n+ 1 v n } . 

Identify members belonging to the same equivalence classes and replace them 
by their canonical representatives. For example, u n v n -\ = u n+ \v„. Adhering to 
lexicographical ordering, use u n v n -\ instead of u n+ \v n . Doing so, replace 5? by 
8? = { , u n v„ _ i , u„ v„ } , which has the building blocks of the density. Linearly com- 
bine the monomials in ST with undetermined coefficients c,- to get the candidate 
density of rank 3 : 

p =c\ ul + c 2 u„v„- 1 +c 3 m„v„. (28) 



3.2 Compute the Undetermined Coefficients Ci 

Compute D t p and use (0) to eliminate u„ and v n and their shifts. Next, introduce the 
main representatives to get 

E = (3c 1 -C2)m 2 v„_i + (c3-3c 1 )m 2 v„ + (c3-c 2 )v„v„ +1 

+ - c 3 )u n u n+l v„ + (c 2 - c 3 )\'l+AJ, (29) 

with 

J = (c 3 - c 2 )v„_iv„ + C 2 M„_lM„V„_ 1 + c 2 \'l_i- (30) 
Set E — AJ = to get the linear system 

3ci — c 2 = 0, C3 — 3ci = 0, c 2 — C3 = 0. (31) 

Select ci = I and substitute the solution c\ = i,c 2 = C3 = 1, into ( f28l and (l30l to 
obtain p' 3 ' in (fT2l with matching flux /' 3 ' = m„_i«„v„_i + v 2 _ t . 
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As an example, we will now compute the symmetry G*- 2 -* = (Gj 2 ' G 2 2 ') T with 
rankG = (3 4) T given in (l20t , 



4.1 Construct the Form of the Symmetry 

Listing all monomials in u„ and v n of ranks 3 and 4, or less: 

^f 2 = {ut,ul,ulv n ,ul,u n v n ,u n ,v 2 n ,v„}. 

Next, for each monomial in «ifi and J^, introduce the necessary f -derivatives so 
that each term exactly has ranks 3 and 4, respectively. At the same time, use (01 to 
remove all f — derivatives. Doing so, based on Jzfi , 

Q C 3\ 3 

dt' 
d° 
dt' 

d , 

— (k„) = 2m„m„ = 2u„v n - 1 -2m„v„, (32) 



^v"/ — v « — u n*n ~ Un+lVrn 

= J t (Un) = ^(V«-1 -V„) 

= M„_lV n _l — U n V n -l — U n V„ + U n+ \V n . 



Put the terms from the right hand sides of (1321 into a set: 

#i = {(^,M„_lV„_ 1 ,M„V„_l,M„V„,M„ + lV„}. 

Similarly, based on the monomials in Jz?2, construct 

5^2 = {u^ u l-l v n-l,U n -lUnVn-l,ulv n -i,V n -2Vn-hvl-l, 
u\v n , M„M„+ 1 V„ , M 2 + ! V„ , V„- 1 V„ , V 2 , V„ V„ + 1 } . 

Linearly combine the monomials in W\ and #2 with undetermined coefficients c,- to 
get the form of the components of the candidate symmetry: 
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c\ ul+c 2 u n -\v n -\ +c 3 m„v„_i +c 4 u n v n +c 5 u n+ \v n , 

(33) 

C6 u\ + ci u 2 n _ ! v„_ i + c 8 u n - 1 u n v„_ ! + c 9 v„_ i 

+ Cl0V n _ 2 V n -l +C 11 V^_ 1 +C 12 M^'„+C 13 M„M„ +1 V„ 
+ C U ul+l v n + C 15V„-1V„ + C l6 V 2 n + C l7 V„V n+l . 

4.2 Compute the Undetermined Coefficients ci 

To determine the coefficients c, , require that ([T5j holds on any solution of ([T). Com- 
pute D,G and use (Q~|) to remove all u,,_i,u n ,u„+i, etc. Compute the Frechet deriva- 
tive dT9b and, in view of ( fT5l ). equate the resulting expressions. Treat as independent 
all the monomials in u„ and their shifts, to obtain the linear system that determines 
the coefficients c,-. 

Apply the strategy to (0 with fl33l l, to get 

Ci = C6 = Cj = eg = eg = Cio = Cn = C13 = C\(, = 0, 
-C2 = — C3 = C 4 = C 5 = -Ci2 = CM = -Ci5 = Ci 7 . 

Set en = 1 and substitute (O into ([33j to get G (2) = {g[ 2) G^) t , as given in d20b . 

To show how our algorithm filters out completely integrable cases among param- 
eterized systems of DDEs, consider 

it„ = a v„_i-v„, 

v„ = v„ (j3 u„ -m„+i), (34) 

where a and j3 are nonzero constant parameters. lfl"9l have shown that (l34l is com- 
pletely integrable if and only if a = J3 = 1 . 

Using our algorithm, one can easily compute the compatibility conditions for a 
and j3 so that (l34l admits a polynomial symmetry, say, of rank (3,4). The steps are 
as above, however, the linear system for the c, is parameterized by a and j3 and must 
be analyzed carefully (with, e.g., Grobner basis methods). This analysis leads to the 
condition a = J3 = 1. Details are given in |9( and iTTOl . 



5 Algorithm for Recursion Operators 

We will now construct the recursion operator (l2ol for d5]. In this case all the terms 
in (fJTJ are 2 x 2 matrix operators. 
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G\ 



(2) 



Gf = 
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5.1 Determine the Rank of the Recursion Operator 



The difference in the ranks of symmetries is used to compute the rank of the ele- 
ments of the recursion operator. Use (|6), (l20t and (l2(Tt to compute 

rankG^^Q, rankG^ = (] V (35) 

Assume that MG^ l > = G^ and use the formula 

rank^ i7 =rankGf +1) -rankG^ , (36) 
to compute a rank matrix associated to the operator £ff : 

rank^= Q i) ■ (37) 



5.2 Determine the Form of the Recursion Operator 

We assume that Si = &q + S£\ , where SIq is a sum of terms involving D 1 , 1, and D. 
(The form of 8&\ will be discussed below.) The coefficients of these terms are ad- 
missible power combinations of u„,u n+ i, v„, and v„_i (which come from the terms 
on the right hand sides of (01), so that all the terms have the correct rank. The maxi- 
mum up-shift and down-shift operator that should be included can be determined by 
comparing two consecutive symmetries. Indeed, if the maximum up-shift in the first 
symmetry is u n+p and the maximum up-shift in the next symmetry is u„ +p+r , then 
the associated piece that goes into S?q must have D,D 2 , . . . ,D r . The same argument 
determines the minimum down-shift operator to be included. For ||3}, get 



(^o)ll (#o)l2 
(^0)21 (^0)22 



(38) 



with 



(^b)ll = (c\U n + C2U n+ l)l, 

(^0)12 -c 3 D- 1 +c 4 I, 
^0)21 

+ (cioul + cnu n u n+ i+c 12 ul +1 +Ci 3 V„_i +C 14 V„)D, 

?o)l2 = (ci 5 M„+C 16 M n+I )I. 



h\ = {c 5 ul + c 6 u„u„ +l +c 7 ul +l +c 8 v„_i + c 9 v„)I (39) 



As shown for the continuous case lfl4l . S%\ is a linear combination (with undeter- 
mined coefficients of all suitable products of symmetries and covariants, i.e., 
Frechet derivatives of densities, sandwiching (D — I) -1 . Hence, 
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^^(D-r 1 ®/, (40) 

j k 

where ® denotes the matrix outer product, defined as 



Gf 



GW(D-I)-VS' G^(D-I)-Vif 
G^D-I)- 1 ^' G W(D-I)-VS' 



(41) 



Only the pair (G^^pi"'') can be used, otherwise the ranks in (|37| | would be ex- 
ceeded. Use ( TBl i and dT9b , to compute 

Pi 0) '=f0 -LlV (42) 



(43) 



From (1401 , after renaming cjo to cn, obtain 

ci 7 (v„_i-v„)(D-I)- 1 il 

ci 7 v„(m„-m„ + i)(D-I) _1 ^1 



Add (O and g3), to get 



with 



= £), (44) 



= (C1M„ +C 2 M„ + l)I, 

*12 = C3D- 1 +c 4 I + ci 7 (v„_i-v„)(D-I)- 1 -I, 

^21 = (c 5 ul + c 6 u n u n+ i+c 7 ul +1 +c s v n -i+c 9 v n )l (45) 

+ (C10M^+C11M„M„+1 +Cl2^ + l+Cl 3 V„_l +Cl 4 V„)D, 

&22 = (ci5M„ + c 16 m„ +1 )I + c 17 v„(m„-m„ +1 )(D-I) _1 — I. 



5.3 Determine the unknown coefficients 



Compute all the terms in d2TT > to find the c,. Refer to lfl5l for the details of the 
computation, resulting in c 2 = C5 = = c-j = cs = cio = cn = C12 = C13 = C15 = 0, 
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and c\ = C3 = C4 = eg = C14 = cig = 1, and c 17 = — 1. Substitute the constants into 
(l44l to get (|26]). 



6 Conclusions and Future Research 

In this paper we presented algorithms for the symbolic computation of polynomial 
conservation laws, generalized symmetries, and recursion operators for systems of 
nonlinear DDEs. We used the Toda lattice to illustrate the steps of the algorithms. 
The algorithms have been implemented in Mathematica and can be used to test the 
complete integrability of nonlinear DDEs. 

Although our algorithm successfully finds conservation laws, generalized sym- 
metries, and recursion operators for various Volterra and Toda lattices as well as the 
Ablowitz-Ladik lattice, the current recursion operator algorithm fails on nonlinear 
DDEs due to Belov and Chaltikian and Blaszak and Marciniak. In future research we 
intend to generalize the recursion operator algorithm so that it can cover a broader 
class of lattices. 

Acknowledgements J.A. Sanders, J. -P. Wang, M. Hickman, and B. Deconinck are gratefully ac- 
knowledged for valuable discussions. 
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